void InitWatAngVels( real extTemperature, VecR ssdInertia, FILE *fPtr ) 
{ /* derived from allen/tildesley p. 170*/
  VecR e, w, angMomentum, totAngMom; 
  int n; real f;
  fprintf( fPtr, "InitWatAngVels... " ); fflush(fPtr);
  VZero( totAngMom );
  //  printf("ssdInertia = %f %f %f\n", ssdInertia.x, ssdInertia.y,ssdInertia.z );
  f = sqrt( 3 * kB_IN_BRAHMS_UNITS * extTemperature / ( VCSum( ssdInertia ) ) );
  //  printf("f = %f\n", f);
  DO_SITE {
    if ( site[ n ].mechType == NSRB ) {
      VCross( angMomentum, siteCoordinates[ n ], siteVelocity[ n ] );
      VVAdd( totAngMom , angMomentum );
      VRand( &e );
      /*       f = 0.5 * velMag / sqrt (VWLenSq (ssdInertia, e)); */
      VSCopy( w, f, e );                         /* body-fixed angular velocity */ 
      VMul( angMomentum, ssdInertia, w );        /* body-fixed angular momentum */
      VCopy( siteAngularMomentum[ n ], angMomentum ); /* body-fixed angular momenta */
      MVMulT( angMomentum, siteRotationMatrix[ n ].u, siteAngularMomentum[ n ] ); /* u^S = Q^T * u^b*/
      VVAdd( totAngMom, angMomentum );
    }
  }
  
  /*   DO_SITE { /\* remove total ang mom *\/ */
  /*       VVSAdd( siteAngularMomentum[ n ], - 1.0 / nSites, totAngMom ); */
  /*     } */

  /*   /\* *** compute total angular momentum; assumption: system COM is the origin *\/ */
  /*   VZero( totAngMom ); */
  /*   DO_SITE { */
  /*     VCross( angMomentum, site[ n ], siteVelocity[ n ] ); */
  /*     VVAdd( totAngMom , angMomentum ); */
  /*     MVMulT( angMomentum, siteRotationMatrix[ n ].u, siteAngularMomentum[ n ] ); /\* u^S = Q^T * u^b*\/ */
  /*     VVAdd( totAngMom, angMomentum ); */
  /*   } */
  fprintf( fPtr, "done\n" ); fflush(fPtr);
}
